#ifndef UHM_SOLVER_OPS_TREE_HXX_
#define UHM_SOLVER_OPS_TREE_HXX_

namespace uhm {

  template<class S_> static Int_ op_uhm_branch_BXR(S_ *s, typename S_::_Tree_::Tnd_* me,
                                                   Int_ b, Int_ x, Int_ r);

  template<class S_> static Int_ op_uhm_branch_BXR(S_ *s, typename S_::_Tree_::Tnd_* me,
                                                   Int_ b, Int_ x, Int_ r) {
    typename S_::_Tree_::_Data_ elt = me->get_data();
    typename S_::_Mesh_::_Uhm_ &uhm = elt->uhm();

    if (me->is_leaf()) {

      // Update to work cache
      Map_< typename S_::_Mesh_::_Id_Nd_> factor, schur;
      elt->nds2map(factor, schur);

      typename S_::_Cache_::Entry_ in(factor, schur), *out;
      UHM_ERROR(s->work().pull(in, &out), UHM_ERROR_SOLVER_PULLED_WORK);

      Mapper_< typename S_::_Mesh_::_Id_Nd_ >
        m_s2f(&out->schur(), &elt->factor()),
        m_s2s(&out->schur(), &elt->schur() );

      m_s2f.relate();
      m_s2s.relate();

      // Mapping B - Do not branch out to B in Work. 

      // Mapping X
      if (x) {
        copy_reverse< typename S_::_Mesh_::_Id_Nd_ >(m_s2f, out->uhm().XB(), uhm.XT());
        copy_reverse< typename S_::_Mesh_::_Id_Nd_ >(m_s2s, out->uhm().XB(), uhm.XB());
      }

      // Mapping R
      if (r) {
        copy_reverse< typename S_::_Mesh_::_Id_Nd_ >(m_s2f, out->uhm().RB(), uhm.RT());
        copy_reverse< typename S_::_Mesh_::_Id_Nd_ >(m_s2s, out->uhm().RB(), uhm.RB());
      }

    } else {
      for (Int_ i=0;i<me->get_n_children();++i) {
        typename S_::_Tree_::_Data_ c_elt = me->get_child(i)->get_data();

        Mapper_< typename S_::_Mesh_::_Id_Nd_ > m_s2f, m_s2s;
        c_elt->relation().get(m_s2f, m_s2s);
        
        typename S_::_Mesh_::_Uhm_ &c_uhm = c_elt->uhm();

        Map_< ID_Node_H_<1> > c_s2f, e_s2f, c_s2s, e_s2s;
        m_s2f.coarse(c_s2f, e_s2f);
        m_s2s.coarse(c_s2s, e_s2s);

        Mapper_< ID_Node_H_<1> > s2f(&c_s2f, &e_s2f), s2s(&c_s2s, &e_s2s);
        s2f.relate();
        s2s.relate();

        // Mapping B - Do not branch out to B.
        
        // Mapping X
        if (x) {
          copy_reverse< ID_Node_H_<1> >(s2f, c_uhm.XB(), uhm.XT());
          copy_reverse< ID_Node_H_<1> >(s2s, c_uhm.XB(), uhm.XB());
        }

        // Mapping R
        if (r) {
          copy_reverse< ID_Node_H_<1> >(s2f, c_uhm.RB(), uhm.RT());
          copy_reverse< ID_Node_H_<1> >(s2s, c_uhm.RB(), uhm.RB());
        }
        
        // After update children, leave the matrices
      }
    }

    return true;
  }

  template<class S_> Int_ op_uhm_branch_X(S_ *s, typename S_::_Tree_::Tnd_* me) {
    return op_uhm_branch_BXR(s, me, 0, 1, 0);
  }

  template<class S_> Int_ op_uhm_branch_R(S_ *s, typename S_::_Tree_::Tnd_* me) {
    return op_uhm_branch_BXR(s, me, 0, 0, 1);
  }

}


#endif
